PEC 3 - Svetlana Bazhenova y Pablo Castejon Molina
2022-06-11
Apartado 1
El set de datos elegido para este trabajo es el “Goldman Osteometric Data Set”. Este dataset incluye datos relativos a 1538 esqueletos humanos de distintas epocas del Holoceno. Los datos recogidos comprenden variables lineales de 4 huesos largos de las piernas (humero, radio, femur y tibia) y tres variables extra de la pelvis. De manera adicional, se estimo el sexo a partir de la morfologia pelvica y dicha estimacion se incluyo en el Data Set.
El Data Set puede descargarse desde la siguiente pagina: https://web.utk.edu/~auerbach/GOLD.htm#gold
Referencias:
Auerbach BM, & Ruff CB. 2004. Human body mass estimation: a comparison of “morphometric” and “mechanical” methods. American Journal of Physical Anthropology 125:331-342.
Auerbach BM, & Ruff CB. 2006. Limb bone bilateral asymmetry: variability and commonality among modern humans. Journal of Human Evolution 50:203-218.
Temple DH, Auerbach BM, Nakatsukasa M, & Larsen CS. 2008. Variation in body proportions between Jomon foragers and Yayoi agriculturalists from prehistoric Japan. American Journal of Physical Anthropology 137:164-174.
Betti L, Von Cramon-Taubadel N, Lycett SJ. 2012. Human pelvis and long bones reveal differential preservation of ancient population history and migration out of Africa. Human Biology 84:139-152.
Plavcan JM. 2012. Body size, size variation, and sexual size dimorphism in Early Homo. Current Anthropology 53:S409-S423.
Stock JT. 2013. The skeletal phenotype of “Negritos” from the Andaman Islands and Philippines relative to global variation among hunter-gatherers. Human Biology 85:67-94.
Roseman CC & Auerbach BM. 2015. Ecogeography, genetics, and the evolution of human body form. Journal of Human Evolution 78:80-90.
Savell KRR, Auerbach BM, & Roseman CC. 2016. Constraint, natural selection, and the evolution of human body form. Proceedings of the National Academy of Sciences USA 113:9492-9497.
Hemos seleccionado estos datos por el interés científico que consideramos que tienen, además de por su gran valor desde el punto de vista de la Bioestádistica: Esta formado por una muestra poblacional muy extensa, incluye múltiples variables de distinto tipo y puede permitirnos inferir mucha información que antes era desconocida a través de un acercamiento y estudio apropiado.
Apartado 2
El archivo se ha descargado de la dirección citada anteriormente, en formato csv, y se nombra como “data”.
download.file("https://web.utk.edu/~auerbach/Goldman.csv", destfile = "data")Se carga el archivo en R y se muestran las primeras lineas:
data1 <- read.csv("data")
head(data1)## Inst ID Sex Age NOTE Location
## 1 AMNH 99.1/69 0 40-50 Tigara - Point Hope, AK Alaska, United States
## 2 AMNH 99.1/70 1 50+ Norton - Point Hope, AK Alaska, United States
## 3 AMNH 99.1/75A 0 30-50 Ipituaq - Point Hope, AK Alaska, United States
## 4 AMNH 99.1/80 0 25-30 Ipituaq - Point Hope, AK Alaska, United States
## 5 AMNH 99.1/83A 1 50+ Ipituaq - Point Hope, AK Alaska, United States
## 6 AMNH 99.1/86B 1 30-40 Ipituaq - Point Hope, AK Alaska, United States
## Element. LHUM RHUM LRAD RRAD LFEM RFEM LTIB RTIB OSCX Metrics. LHML LHEB
## 1 NA 0 1 0 1 0 0 1 0 0 NA 308.5 64
## 2 NA 1 1 1 1 0 0 0 0 0 NA NA NA
## 3 NA 0 0 1 0 0 0 0 0 1 NA 311.0 59
## 4 NA 0 0 0 0 0 0 0 0 0 NA 289.0 55
## 5 NA 0 0 0 0 0 0 0 0 0 NA 295.0 54
## 6 NA 0 0 0 1 1 0 0 0 0 NA 270.5 55
## LHHD LHMLD LHAPD RHML RHEB RHHD RHMLD RHAPD LRML LRMLD LRAPD RRML RRMLD
## 1 47.55 24.30 22.00 NA NA NA NA NA 229.0 14.84 12.02 NA NA
## 2 NA NA NA NA NA NA NA NA NA NA NA NA NA
## 3 44.44 22.20 22.12 310 59.5 44.11 21.20 22.68 NA NA NA 221.5 17.26
## 4 42.94 20.83 20.36 298 58.5 44.41 21.36 22.09 223.5 13.89 11.12 227.0 14.16
## 5 42.51 19.34 19.35 302 55.0 42.06 20.21 19.97 208.0 14.38 9.62 205.0 15.63
## 6 39.74 19.22 19.42 281 55.5 39.84 19.74 19.38 196.5 13.66 8.85 NA NA
## RRAPD LFML LFBL LFEB LFAB LFHD LFMLD LFAPD RFML RFBL RFEB RFAB RFHD
## 1 NA 443 441.5 82.0 70.68 44.92 29.77 31.39 452.0 451.0 85.0 75.85 48.36
## 2 NA 386 383.0 72.0 61.34 44.64 24.97 23.86 383.0 380.5 73.0 62.00 44.09
## 3 11.03 415 410.0 75.0 67.34 43.95 27.64 28.85 415.5 410.0 77.5 67.35 44.17
## 4 11.22 398 392.0 78.5 66.71 43.97 24.82 30.73 400.0 397.0 77.5 66.95 43.77
## 5 9.95 395 393.0 69.0 57.72 40.58 24.69 28.26 396.0 393.0 70.0 57.88 41.13
## 6 NA NA NA NA NA NA NA NA 375.0 372.5 68.0 58.90 40.80
## RFMLD RFAPD LTML LTPB LTMLD LTAPD RTML RTPB RTMLD RTAPD BIB LIBL RIBL
## 1 28.94 33.01 NA NA NA NA 365.5 77.0 22.41 26.75 274.0 162 161
## 2 24.11 24.41 283 NA 19.12 23.73 281.5 NA 20.04 21.35 266.0 NA NA
## 3 27.18 28.39 330 73 21.72 26.96 331.5 73.0 21.07 25.36 NA NA NA
## 4 25.18 31.50 312 69 21.76 24.56 320.5 71.5 21.98 24.41 240.5 NA 145
## 5 24.79 26.88 306 64 19.44 23.43 305.0 NA 20.30 21.45 285.5 160 155
## 6 23.76 24.40 298 65 17.75 24.64 293.5 62.5 17.05 24.15 NA NA NA
## LAcH RAcH Derived. Brachial Crural IL.UL.LL IL.LL.UL CBR.FHD McH.FHD
## 1 52.56 52.24 NA 0.7423015 0.8167598 0.6611316 1.512558 65.64622 64.52696
## 2 49.21 50.18 NA NA 0.7340702 NA NA 65.27654 59.43324
## 3 NA NA NA 0.7133655 0.7965081 0.7131367 1.402256 59.28161 58.75034
## 4 50.30 49.80 NA 0.7674617 0.7926065 0.7252709 1.378795 58.81290 58.32493
## 5 45.60 44.83 NA 0.6917923 0.7724399 0.7203994 1.388119 57.61281 51.57435
## 6 45.70 47.04 NA 0.7126020 0.7886667 0.7040626 1.420328 57.49272 51.45120
## GRINE.FHD AVG.FHD
## 1 69.27952 66.48423
## 2 64.11982 62.94320
## 3 63.42808 60.48668
## 4 62.99716 60.04500
## 5 56.15914 55.11543
## 6 56.03440 54.99277Tras ello, se muestra un resumen de las variables y su tipo:
summary(data1)## Inst ID Sex Age
## Length:1538 Length:1538 Length:1538 Length:1538
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## NOTE Location Element. LHUM
## Length:1538 Length:1538 Mode:logical Min. :0.0000
## Class :character Class :character NA's:1538 1st Qu.:0.0000
## Mode :character Mode :character Median :0.0000
## Mean :0.1008
## 3rd Qu.:0.0000
## Max. :1.0000
##
## RHUM LRAD RRAD LFEM
## Min. :0.00000 Min. :0.0000 Min. :0.0000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.00000
## Median :0.00000 Median :0.0000 Median :0.0000 Median :0.00000
## Mean :0.08127 Mean :0.1404 Mean :0.1268 Mean :0.07022
## 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.00000
## Max. :1.00000 Max. :1.0000 Max. :1.0000 Max. :1.00000
##
## RFEM LTIB RTIB OSCX
## Min. :0.00000 Min. :0.00000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.00000 Median :0.00000 Median :0.00000 Median :0.00000
## Mean :0.06437 Mean :0.08648 Mean :0.08648 Mean :0.02016
## 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.00000
## Max. :1.00000 Max. :1.00000 Max. :1.00000 Max. :1.00000
##
## Metrics. LHML LHEB LHHD LHMLD
## Mode:logical Min. :229.5 Min. :36.51 Min. :29.58 Min. :12.01
## NA's:1538 1st Qu.:287.4 1st Qu.:53.00 1st Qu.:39.76 1st Qu.:18.23
## Median :303.5 Median :57.50 Median :42.72 Median :19.86
## Mean :303.8 Mean :57.44 Mean :42.74 Mean :19.88
## 3rd Qu.:319.0 3rd Qu.:61.00 3rd Qu.:45.72 3rd Qu.:21.46
## Max. :376.0 Max. :75.00 Max. :56.33 Max. :27.22
## NA's :162 NA's :184 NA's :170 NA's :162
## LHAPD RHML RHEB RHHD
## Min. :13.00 Min. :234.0 Min. :42.00 Min. :32.24
## 1st Qu.:18.03 1st Qu.:291.0 1st Qu.:54.00 1st Qu.:39.90
## Median :19.64 Median :307.0 Median :58.00 Median :43.02
## Mean :19.72 Mean :307.6 Mean :58.19 Mean :43.00
## 3rd Qu.:21.30 3rd Qu.:323.5 3rd Qu.:62.00 3rd Qu.:46.00
## Max. :27.29 Max. :383.0 Max. :75.00 Max. :55.67
## NA's :162 NA's :135 NA's :154 NA's :142
## RHMLD RHAPD LRML LRMLD
## Min. :12.25 Min. :13.84 Min. :179.0 Min. : 8.74
## 1st Qu.:18.69 1st Qu.:18.73 1st Qu.:219.0 1st Qu.:12.78
## Median :20.30 Median :20.46 Median :233.0 Median :14.05
## Mean :20.36 Mean :20.46 Mean :233.1 Mean :14.10
## 3rd Qu.:22.04 3rd Qu.:22.19 3rd Qu.:247.0 3rd Qu.:15.32
## Max. :30.44 Max. :27.19 Max. :290.5 Max. :22.15
## NA's :135 NA's :136 NA's :217 NA's :217
## LRAPD RRML RRMLD RRAPD LFML
## Min. : 7.37 Min. :180 Min. : 9.16 Min. : 7.88 Min. :345.0
## 1st Qu.:10.23 1st Qu.:221 1st Qu.:13.12 1st Qu.:10.33 1st Qu.:404.5
## Median :11.11 Median :235 Median :14.40 Median :11.27 Median :428.0
## Mean :11.19 Mean :235 Mean :14.49 Mean :11.34 Mean :427.1
## 3rd Qu.:12.12 3rd Qu.:248 3rd Qu.:15.74 3rd Qu.:12.31 3rd Qu.:448.5
## Max. :16.28 Max. :291 Max. :23.22 Max. :15.68 Max. :531.0
## NA's :217 NA's :201 NA's :198 NA's :198 NA's :117
## LFBL LFEB LFAB LFHD
## Min. :309.5 Min. :58.00 Min. :49.99 Min. :33.64
## 1st Qu.:401.0 1st Qu.:72.00 1st Qu.:62.01 1st Qu.:40.61
## Median :424.0 Median :76.00 Median :66.61 Median :43.54
## Mean :423.5 Mean :76.07 Mean :66.36 Mean :43.43
## 3rd Qu.:445.0 3rd Qu.:81.00 3rd Qu.:70.82 3rd Qu.:46.15
## Max. :530.0 Max. :93.50 Max. :83.68 Max. :57.39
## NA's :122 NA's :158 NA's :160 NA's :117
## LFMLD LFAPD RFML RFBL
## Min. :17.90 Min. :18.88 Min. :341.5 Min. :279.0
## 1st Qu.:23.83 1st Qu.:25.12 1st Qu.:403.0 1st Qu.:399.0
## Median :25.73 Median :27.16 Median :426.0 Median :422.0
## Mean :25.77 Mean :27.17 Mean :425.7 Mean :421.6
## 3rd Qu.:27.65 3rd Qu.:29.18 3rd Qu.:447.5 3rd Qu.:443.5
## Max. :34.54 Max. :37.06 Max. :532.5 Max. :531.0
## NA's :115 NA's :115 NA's :112 NA's :115
## RFEB RFAB RFHD RFMLD
## Min. :58.00 Min. :49.83 Min. :32.97 Min. :17.97
## 1st Qu.:72.00 1st Qu.:61.97 1st Qu.:40.51 1st Qu.:23.54
## Median :76.00 Median :66.19 Median :43.52 Median :25.37
## Mean :76.26 Mean :66.28 Mean :43.50 Mean :25.44
## 3rd Qu.:81.00 3rd Qu.:70.61 3rd Qu.:46.40 3rd Qu.:27.24
## Max. :94.00 Max. :82.85 Max. :57.86 Max. :33.84
## NA's :152 NA's :148 NA's :103 NA's :114
## RFAPD LTML LTPB LTMLD
## Min. :18.34 Min. :276.0 Min. :52.00 Min. :15.22
## 1st Qu.:25.10 1st Qu.:333.0 1st Qu.:65.00 1st Qu.:19.78
## Median :27.14 Median :353.0 Median :69.50 Median :21.37
## Mean :27.24 Mean :353.1 Mean :69.34 Mean :21.46
## 3rd Qu.:29.29 3rd Qu.:372.0 3rd Qu.:74.00 3rd Qu.:23.03
## Max. :37.53 Max. :446.0 Max. :86.00 Max. :31.00
## NA's :114 NA's :135 NA's :186 NA's :139
## LTAPD RTML RTPB RTMLD
## Min. :19.00 Min. :237.0 Min. :50.00 Min. :15.15
## 1st Qu.:24.18 1st Qu.:332.0 1st Qu.:65.00 1st Qu.:20.25
## Median :26.23 Median :353.0 Median :69.00 Median :21.89
## Mean :26.35 Mean :352.4 Mean :69.34 Mean :21.99
## 3rd Qu.:28.54 3rd Qu.:372.0 3rd Qu.:74.00 3rd Qu.:23.79
## Max. :37.94 Max. :444.0 Max. :85.00 Max. :29.83
## NA's :140 NA's :138 NA's :189 NA's :143
## RTAPD BIB LIBL RIBL LAcH
## Min. :18.59 Min. :184.0 Min. :105 Min. :106.0 Min. :38.64
## 1st Qu.:23.38 1st Qu.:251.0 1st Qu.:145 1st Qu.:145.0 1st Qu.:45.99
## Median :25.21 Median :263.0 Median :151 Median :151.0 Median :49.05
## Mean :25.37 Mean :262.4 Mean :151 Mean :151.1 Mean :48.87
## 3rd Qu.:27.32 3rd Qu.:274.0 3rd Qu.:158 3rd Qu.:158.0 3rd Qu.:51.72
## Max. :35.20 Max. :324.0 Max. :181 Max. :189.0 Max. :62.94
## NA's :144 NA's :69 NA's :359 NA's :361 NA's :167
## RAcH Derived. Brachial Crural
## Min. :37.89 Mode:logical Min. :0.6779 Min. :0.6920
## 1st Qu.:46.05 NA's:1538 1st Qu.:0.7446 1st Qu.:0.8092
## Median :49.03 Median :0.7653 Median :0.8284
## Mean :48.92 Mean :0.7653 Mean :0.8277
## 3rd Qu.:51.80 3rd Qu.:0.7857 3rd Qu.:0.8459
## Max. :62.99 Max. :1.0769 Max. :0.9654
## NA's :163 NA's :75 NA's :48
## IL.UL.LL IL.LL.UL CBR.FHD McH.FHD
## Min. :0.5993 Min. :1.299 Min. :36.09 Min. :34.69
## 1st Qu.:0.6803 1st Qu.:1.421 1st Qu.:54.61 1st Qu.:50.99
## Median :0.6916 Median :1.446 Median :60.05 Median :57.47
## Mean :0.6922 Mean :1.446 Mean :60.15 Mean :57.41
## 3rd Qu.:0.7036 3rd Qu.:1.470 3rd Qu.:65.20 3rd Qu.:63.74
## Max. :0.7696 Max. :1.669 Max. :92.75 Max. :89.12
## NA's :120 NA's :120 NA's :19 NA's :19
## GRINE.FHD AVG.FHD
## Min. :39.06 Min. :38.30
## 1st Qu.:55.57 1st Qu.:53.74
## Median :62.14 Median :59.83
## Mean :62.07 Mean :59.88
## 3rd Qu.:68.48 3rd Qu.:65.68
## Max. :94.19 Max. :92.02
## NA's :19 NA's :19str(data1)## 'data.frame': 1538 obs. of 69 variables:
## $ Inst : chr "AMNH" "AMNH" "AMNH" "AMNH" ...
## $ ID : chr "99.1/69" "99.1/70" "99.1/75A" "99.1/80" ...
## $ Sex : chr "0" "1" "0" "0" ...
## $ Age : chr "40-50" "50+" "30-50" "25-30" ...
## $ NOTE : chr "Tigara - Point Hope, AK" "Norton - Point Hope, AK" "Ipituaq - Point Hope, AK" "Ipituaq - Point Hope, AK" ...
## $ Location : chr "Alaska, United States" "Alaska, United States" "Alaska, United States" "Alaska, United States" ...
## $ Element. : logi NA NA NA NA NA NA ...
## $ LHUM : int 0 1 0 0 0 0 0 0 0 0 ...
## $ RHUM : int 1 1 0 0 0 0 0 0 0 0 ...
## $ LRAD : int 0 1 1 0 0 0 0 0 0 0 ...
## $ RRAD : int 1 1 0 0 0 1 0 0 0 0 ...
## $ LFEM : int 0 0 0 0 0 1 0 0 0 0 ...
## $ RFEM : int 0 0 0 0 0 0 0 0 0 0 ...
## $ LTIB : int 1 0 0 0 0 0 0 0 0 0 ...
## $ RTIB : int 0 0 0 0 0 0 0 0 0 0 ...
## $ OSCX : int 0 0 1 0 0 0 1 0 0 0 ...
## $ Metrics. : logi NA NA NA NA NA NA ...
## $ LHML : num 308 NA 311 289 295 ...
## $ LHEB : num 64 NA 59 55 54 55 63 NA 63.5 57 ...
## $ LHHD : num 47.5 NA 44.4 42.9 42.5 ...
## $ LHMLD : num 24.3 NA 22.2 20.8 19.3 ...
## $ LHAPD : num 22 NA 22.1 20.4 19.4 ...
## $ RHML : num NA NA 310 298 302 ...
## $ RHEB : num NA NA 59.5 58.5 55 55.5 60 57 62 58 ...
## $ RHHD : num NA NA 44.1 44.4 42.1 ...
## $ RHMLD : num NA NA 21.2 21.4 20.2 ...
## $ RHAPD : num NA NA 22.7 22.1 20 ...
## $ LRML : num 229 NA NA 224 208 ...
## $ LRMLD : num 14.8 NA NA 13.9 14.4 ...
## $ LRAPD : num 12.02 NA NA 11.12 9.62 ...
## $ RRML : num NA NA 222 227 205 ...
## $ RRMLD : num NA NA 17.3 14.2 15.6 ...
## $ RRAPD : num NA NA 11.03 11.22 9.95 ...
## $ LFML : num 443 386 415 398 395 ...
## $ LFBL : num 442 383 410 392 393 ...
## $ LFEB : num 82 72 75 78.5 69 NA 79 73 79 73 ...
## $ LFAB : num 70.7 61.3 67.3 66.7 57.7 ...
## $ LFHD : num 44.9 44.6 44 44 40.6 ...
## $ LFMLD : num 29.8 25 27.6 24.8 24.7 ...
## $ LFAPD : num 31.4 23.9 28.9 30.7 28.3 ...
## $ RFML : num 452 383 416 400 396 ...
## $ RFBL : num 451 380 410 397 393 ...
## $ RFEB : num 85 73 77.5 77.5 70 68 80 74 81 71 ...
## $ RFAB : num 75.8 62 67.3 67 57.9 ...
## $ RFHD : num 48.4 44.1 44.2 43.8 41.1 ...
## $ RFMLD : num 28.9 24.1 27.2 25.2 24.8 ...
## $ RFAPD : num 33 24.4 28.4 31.5 26.9 ...
## $ LTML : num NA 283 330 312 306 298 351 361 338 311 ...
## $ LTPB : num NA NA 73 69 64 65 73 69 74.5 68 ...
## $ LTMLD : num NA 19.1 21.7 21.8 19.4 ...
## $ LTAPD : num NA 23.7 27 24.6 23.4 ...
## $ RTML : num 366 282 332 320 305 ...
## $ RTPB : num 77 NA 73 71.5 NA 62.5 73 70 71 70 ...
## $ RTMLD : num 22.4 20 21.1 22 20.3 ...
## $ RTAPD : num 26.8 21.4 25.4 24.4 21.4 ...
## $ BIB : num 274 266 NA 240 286 ...
## $ LIBL : num 162 NA NA NA 160 NA NA 158 160 160 ...
## $ RIBL : num 161 NA NA 145 155 NA NA 156 161 160 ...
## $ LAcH : num 52.6 49.2 NA 50.3 45.6 ...
## $ RAcH : num 52.2 50.2 NA 49.8 44.8 ...
## $ Derived. : logi NA NA NA NA NA NA ...
## $ Brachial : num 0.742 NA 0.713 0.767 0.692 ...
## $ Crural : num 0.817 0.734 0.797 0.793 0.772 ...
## $ IL.UL.LL : num 0.661 NA 0.713 0.725 0.72 ...
## $ IL.LL.UL : num 1.51 NA 1.4 1.38 1.39 ...
## $ CBR.FHD : num 65.6 65.3 59.3 58.8 57.6 ...
## $ McH.FHD : num 64.5 59.4 58.8 58.3 51.6 ...
## $ GRINE.FHD: num 69.3 64.1 63.4 63 56.2 ...
## $ AVG.FHD : num 66.5 62.9 60.5 60 55.1 ...Como puede verse, el dataset tiene una gran cantidad de variables de distintos tipos (1538 observaciones de 69 variables), incluyendo numéricas (las relativas a las mediciones en los huesos), cualitativas (las relativas a la presencia o ausencia de uno de los huesos) y de caracteres (información adicional como el sexo, el código de identificación o el lugar del yacimiento).
La descripcion de las medidas se encuentra en el PDF encontrado en el mismo sitio que los datos: https://web.utk.edu/~auerbach/GoldMeasures.pdf
Podemos ver que las columnas Elements, Metrics y Derived estan compuestas unicamente de NA y sirven unicamente a separar las variables por tipo establecido por los autores (Elements, Metrics y Derived) dentro del dataset.
Otras cosas importantes de notar:
- las variables de tipo Elements (LHUM, RHUM, LRAD, RRAD, LFEM, RFEM, LTIB, RTIB, OSCX) son de tipo categorico (0 y 1) y tienen un significado contraintuitivo (0 significa que el elemento es presente y almenos una medida ha sido prelevada y 1 que es ausente)
- Las variables de tipo Metrics descritas en el PDF son presentes dos veces: una para la izquierda, el nombre de la variable empezando con la letra L (left), una para la derecha, empezando con R (right), por ejemplo en vez de tener senciallamente IBL (descripcion: Maximum Iliac Blade Length), en el dataset tenemos a RIBL y LIBL
Apartado 3
Preguntas posibles:
- ¿Que huesos se preservan mas de la derecha si se preserva la izquierda? (Elements)
- ¿Los museos estan especializados en origen de huesos?
- (Metrics) la variable BIB (Bi‐iliac Breadth) no tiene derecha ni izquierda. Hay alguna variable con cuya ella esta correlada (que no sea genero)
- Variable Sex : podemos hacer una predicion de genero a partir de este dataset para las observaciones donde Sex esta marcado como 0? y 1? ? Los resultados coinciden con el metodo que usaron los autores?
- ¿Tienen de media los húmeros derechos e izquierdos (LHML y RHML) la misma longitud máxima o hay diferencias?
- ¿Que porcentaje de la población total corresponde a los individuos de +50 años? ¿En qué llacimiento hay mas individuos de este rango de población?
Respuetas:
a) ¿Que huesos se preservan mas de la derecha si se preserva la izquierda? (Elements)
Tenemos que mirar la proporcion de preservacion de las variables de tipo Elements (LHUM, RHUM, LRAD, RRAD, LFEM, RFEM, LTIB, RTIB). La poblacion completa para medir el porcentaje de preservacion de los huesos derechas seria para cada par (LHUM y RHUM, LRAD y RRAD, LFEM y RFEM, LTIB y RTIB) la poblacion donde la preservacion de huesos izquierdos seria completa (== 0)
#Se prepara un nuevo data set, donde esta preservado el LHUM
LHUM_preservado <- data1[data1$LHUM == '0', ]
LHUM_preservado <- LHUM_preservado[!is.na(LHUM_preservado$LHUM),]
# Y aqui preparamos el subdataset de LHUM que solo tiene los RHUM preservados
RHUM_preservado <- LHUM_preservado[LHUM_preservado$RHUM == '0', ]
RHUM_preservado <- RHUM_preservado[!is.na(RHUM_preservado$RHUM),]
# proporcion de RHUM preservados cuando LHUM esta preservado
print("humerus: ")## [1] "humerus: "nrow(RHUM_preservado)/nrow(LHUM_preservado)## [1] 0.934201#Se prepara un nuevo data set, donde esta preservado el LRAD
LRAD_preservado <- data1[data1$LRAD == '0', ]
LRAD_preservado <- LRAD_preservado[!is.na(LRAD_preservado$LRAD),]
# Y aqui preparamos el subdataset de LRAD que solo tiene los RRAD preservados
RRAD_preservado <- LRAD_preservado[LRAD_preservado$RRAD == '0', ]
RRAD_preservado <- RRAD_preservado[!is.na(RRAD_preservado$RRAD),]
# proporcion de RRAD preservados cuando LRAD esta preservado
print("radius: ")## [1] "radius: "nrow(RRAD_preservado)/nrow(LRAD_preservado)## [1] 0.8933434#Se prepara un nuevo data set, donde esta preservado el LFEM
LFEM_preservado <- data1[data1$LFEM == '0', ]
LFEM_preservado <- LFEM_preservado[!is.na(LFEM_preservado$LFEM),]
# Y aqui preparamos el subdataset de LFEM que solo tiene los RFEM preservados
RFEM_preservado <- LFEM_preservado[LFEM_preservado$RFEM == '0', ]
RFEM_preservado <- RFEM_preservado[!is.na(RFEM_preservado$RFEM),]
# proporcion de RFEM preservados cuando LFEM esta preservado
print("femur: ")## [1] "femur: "nrow(RFEM_preservado)/nrow(LFEM_preservado)## [1] 0.9426573#Se prepara un nuevo data set, donde esta preservado el LTIB
LTIB_preservado <- data1[data1$LTIB == '0', ]
LTIB_preservado <- LTIB_preservado[!is.na(LTIB_preservado$LTIB),]
# Y aqui preparamos el subdataset de LTIB que solo tiene los RTIB preservados
RTIB_preservado <- LTIB_preservado[LTIB_preservado$RTIB == '0', ]
RTIB_preservado <- RTIB_preservado[!is.na(RTIB_preservado$RTIB),]
# proporcion de RTIB preservados cuando LTIB esta preservado
print("tibia: ")## [1] "tibia: "nrow(RTIB_preservado)/nrow(LTIB_preservado)## [1] 0.9338078Vemos que la proporcion mas alta la tienen LFEM y RFEM que son los femures. No deberia ser una sorpresa ya que es el hueso mas largo y fuerte del cuerpo humano asi que lo logico seria que es el que se preserve mas, también proporcialmente entre la derecha y la izquierda.
b) ¿Los museos estan especializados en origen de huesos?
Podemos hacer unos barplots por museo
par(mar=c(5,10,4,4))
barplot(table(data1$Location,data1$Inst),horiz=TRUE,las=2, cex.names=.6,col = rainbow(length(unique(data1$Location))))
El barplot sale muy colorado. Podemos ver que WOAC, SfAP, KU, KSU y CMNH estan muy especializados cada uno en una Location (la barra es mayoritarmente de un unico color). Cual? Para averiguarlo, podemos hacer un subset del dataset componiendolo unicamente de observaciones de estos museos.
inst_woac <- data1[data1$Inst == 'WOAC', ]
inst_sfap <- data1[data1$Inst == 'SfAP', ]
inst_ku <- data1[data1$Inst == 'KU', ]
inst_ksu <- data1[data1$Inst == 'KSU', ]
inst_cmnh <- data1[data1$Inst == 'CMNH', ]
# mostramos la reparticion de los origenes por museo
table(inst_woac$Location,inst_woac$Inst)##
## WOAC
## Kentucky, United States 61table(inst_sfap$Location,inst_sfap$Inst)##
## SfAP
## Germany 108
## Tanzania 1table(inst_ku$Location,inst_ku$Inst)##
## KU
## Japan 66table(inst_ksu$Location,inst_ksu$Inst)##
## KSU
## Ohio, United States 21table(inst_cmnh$Location,inst_cmnh$Inst)##
## CMNH
## Colorado, United States 3
## Hawaii, United States 2
## Ohio, United States 170Vemos que los museos “especializados” en un origen de los huesos tienden a tener colleciones regionales de la region o pais donde esta el museo. Asi mismo, WOAC – Webb Osteology and Archaeology Collection (que se encuentra en el estado de Kentucky, EEUU) tiene unicamente los huesos con origen en Kentucky. SfAP – Staatssammlung fur Anthropologie und Palaeoanatomie (que es una collecion alemana) los huesos con origen en Alemania (a l’exepcion de un ejemplar de Tanzania), KU – Kyoto University (Kyodai) con origen en Japon, KSU – Kent State University y CMNH – Cleveland Museum of Natural History los huesos con origen de Ohio (que es el estado donde estan).
c) (Metrics) la variable BIB (Bi‐iliac Breadth) no tiene derecha ni izquierda. Hay alguna variable con cuya ella esta correlada (que no sea genero)
#Se importa zoo para realizar la sustitución rápidamente
#install.packages("zoo")
# ponemos las metricas en un dataset a parte
metrics <- data1[,c(18:60)]
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric# reemplazamos los NA en cada variable con la media de la variable
metrics2 <- na.aggregate(metrics)
# miramos la correlacion de cada variable con BIB
cor(metrics2[-39], metrics2$BIB)## [,1]
## LHML 0.4711816
## LHEB 0.4696301
## LHHD 0.4859928
## LHMLD 0.4124235
## LHAPD 0.3888506
## RHML 0.4915438
## RHEB 0.4690396
## RHHD 0.5104892
## RHMLD 0.4360081
## RHAPD 0.4171695
## LRML 0.3725394
## LRMLD 0.3487594
## LRAPD 0.3437768
## RRML 0.3657437
## RRMLD 0.3618544
## RRAPD 0.3654109
## LFML 0.4557349
## LFBL 0.4552652
## LFEB 0.5052740
## LFAB 0.4620532
## LFHD 0.5474541
## LFMLD 0.5451732
## LFAPD 0.4323758
## RFML 0.4606696
## RFBL 0.4523228
## RFEB 0.4911563
## RFAB 0.4374224
## RFHD 0.5463900
## RFMLD 0.5285022
## RFAPD 0.4349595
## LTML 0.3842317
## LTPB 0.4499900
## LTMLD 0.3859132
## LTAPD 0.3545380
## RTML 0.3783401
## RTPB 0.4636810
## RTMLD 0.3782224
## RTAPD 0.3675581
## LIBL 0.6501100
## RIBL 0.6679058
## LAcH 0.5372142
## RAcH 0.5320178La correlacion mas grande de BIB es con RIBL (Right Maximum Iliac Blade Length) (0.6679058) y LIBL (Left Maximum Iliac Blade Length) (0.6501100).
d) Variable Sex : podemos hacer una predicion de genero a partir de este dataset para las observaciones donde Sex esta marcado como 0? y 1? ? Los resultados coinciden con el metodo que usaron los autores?
Podemos usar un algoritmo de clasificacion, por ejempo Support Vector Machine (SVM)
Para aplicar ese algoritmo, hemos seguido en parte las instrucciones encontradas aqui: https://www.geeksforgeeks.org/classifying-data-using-support-vector-machinessvms-in-r/
# separamos el dataset data2 en grupos test y train. El test va a incluir todas las observaciones donde el genero no esta seguro (marcado 0? o 1?) y el train todas las demas observaciones.
data2_test <- data1[data1$Sex == '0?'| data1$Sex == '1?', ]
data2_train <- data1[data1$Sex == '0'| data1$Sex == '1', ]
# en el dataset test reemplazamos los valores 0? y 1? con un el equivalente 0 y 1 para compararlos después de entrenar el algoritmo
data2_test$Sex<-replace(data2_test$Sex, data2_test$Sex=="0?","0")
data2_test$Sex<-replace(data2_test$Sex, data2_test$Sex=="1?","1")
# para el dataset final de train y test, necesitamos unicamente las columnas de metricas
metricsgen_train <- data2_train[,c(18:60)]
metricsgen_test <- data2_test[,c(18:60)]
# para la variable que va a ser predecida, tomamos a parte la columna Sex en ambos dataset train y test
sex_train <- data2_train[,c(3)]
sex_test <- data2_test[,c(3)]
# pasamos a tipo factor la variable predecida en ambos train y test
sex_train<-factor(sex_train, levels = c(0, 1))
sex_test<-factor(sex_test, levels = c(0, 1))
# reemplazamos con la media de la variable las observaciones que tienen NA en ambos train y test
library(zoo)
metricsgen_train <- na.aggregate(metricsgen_train)
metricsgen_test <- na.aggregate(metricsgen_test)
head(metricsgen_test)## LHML LHEB LHHD LHMLD LHAPD RHML RHEB RHHD RHMLD RHAPD LRML LRMLD
## 166 294.0 53 39.17 17.85 18.89 296.0 55 39.48 17.98 19.05 227.5000 13.7800
## 169 290.0 53 37.43 17.85 20.62 294.0 55 36.55 17.65 20.51 226.5000 14.7000
## 171 291.0 56 40.88 20.05 18.71 291.5 59 41.01 21.55 19.61 233.9375 13.8475
## 509 324.5 57 42.05 20.81 18.95 329.0 57 42.15 20.16 18.69 259.0000 15.2400
## 744 293.0 58 36.09 21.05 21.92 300.0 55 40.77 21.34 22.15 233.9375 13.8475
## 994 287.0 50 36.55 16.85 16.64 293.0 52 38.25 18.35 17.25 220.5000 12.6500
## LRAPD RRML RRMLD RRAPD LFML LFBL LFEB LFAB LFHD
## 166 10.4600 236.6875 14.3475 10.83625 416.5 413.0000 72.50000 61.91000 40.75000
## 169 11.1600 230.0000 13.8500 11.38000 411.0 408.5000 69.50000 60.72000 38.54000
## 171 10.7125 230.5000 14.9900 11.04000 403.0 402.0000 71.00000 61.00000 40.46000
## 509 11.6600 263.0000 15.0100 11.31000 460.5 459.0000 80.50000 67.11000 43.27000
## 744 10.7125 236.6875 14.3475 10.83625 424.5 420.6667 73.22222 63.47889 41.23111
## 994 10.2300 228.0000 14.6800 10.74000 411.0 404.5000 72.00000 62.34000 40.44000
## LFMLD LFAPD RFML RFBL RFEB RFAB RFHD RFMLD
## 166 24.33000 27.17000 415.9375 412.1875 73.05556 63.83111 40.4575 24.62333
## 169 23.36000 27.77000 410.0000 407.0000 70.00000 60.17000 38.9100 22.57000
## 171 25.56000 26.64000 404.5000 402.0000 73.00000 61.83000 40.7000 24.65000
## 509 25.38000 27.53000 415.9375 412.1875 79.50000 68.73000 40.4575 27.05000
## 744 24.44111 26.79444 402.0000 398.0000 69.00000 60.93000 36.6900 25.98000
## 994 26.94000 29.01000 403.0000 395.0000 70.00000 60.87000 39.7600 25.56000
## RFAPD LTML LTPB LTMLD LTAPD RTML RTPB RTMLD RTAPD BIB
## 166 27.08889 350.5 64.0 20.62 25.65 357.1111 68.38889 21.41556 24.91889 265.0
## 169 28.18000 359.0 64.5 20.65 25.56 360.5000 65.00000 21.50000 25.39000 266.0
## 171 27.13000 334.0 65.0 22.58 25.67 337.5000 66.00000 22.39000 24.79000 269.5
## 509 26.69000 389.5 76.0 22.00 29.82 391.5000 76.00000 21.40000 29.70000 273.0
## 744 28.19000 371.0 74.0 21.16 24.00 376.0000 68.00000 20.73000 22.79000 257.0
## 994 27.43000 343.0 67.0 19.49 26.14 341.0000 67.50000 20.30000 24.59000 280.5
## LIBL RIBL LAcH RAcH
## 166 140.0000 147.375 46.56000 46.63
## 169 146.0000 145.000 44.44000 45.69
## 171 143.0000 146.000 47.70000 47.12
## 509 150.0000 151.000 46.79111 49.98
## 744 153.0000 154.000 46.06000 44.87
## 994 145.6667 147.375 46.63000 46.44head(metricsgen_train)## LHML LHEB LHHD LHMLD LHAPD RHML RHEB RHHD
## 1 308.5000 64.0000 47.55000 24.30000 22.00000 307.5837 58.21747 43.0266
## 2 303.7822 57.4669 42.76761 19.89165 19.72556 307.5837 58.21747 43.0266
## 3 311.0000 59.0000 44.44000 22.20000 22.12000 310.0000 59.50000 44.1100
## 4 289.0000 55.0000 42.94000 20.83000 20.36000 298.0000 58.50000 44.4100
## 5 295.0000 54.0000 42.51000 19.34000 19.35000 302.0000 55.00000 42.0600
## 6 270.5000 55.0000 39.74000 19.22000 19.42000 281.0000 55.50000 39.8400
## RHMLD RHAPD LRML LRMLD LRAPD RRML RRMLD RRAPD
## 1 20.36844 20.47182 229.0000 14.84000 12.02000 234.9537 14.49474 11.33899
## 2 20.36844 20.47182 233.0636 14.09854 11.19121 234.9537 14.49474 11.33899
## 3 21.20000 22.68000 233.0636 14.09854 11.19121 221.5000 17.26000 11.03000
## 4 21.36000 22.09000 223.5000 13.89000 11.12000 227.0000 14.16000 11.22000
## 5 20.21000 19.97000 208.0000 14.38000 9.62000 205.0000 15.63000 9.95000
## 6 19.74000 19.38000 196.5000 13.66000 8.85000 234.9537 14.49474 11.33899
## LFML LFBL LFEB LFAB LFHD LFMLD LFAPD RFML RFBL
## 1 443.0000 441.5000 82.00000 70.68000 44.92000 29.77000 31.39000 452.0 451.0
## 2 386.0000 383.0000 72.00000 61.34000 44.64000 24.97000 23.86000 383.0 380.5
## 3 415.0000 410.0000 75.00000 67.34000 43.95000 27.64000 28.85000 415.5 410.0
## 4 398.0000 392.0000 78.50000 66.71000 43.97000 24.82000 30.73000 400.0 397.0
## 5 395.0000 393.0000 69.00000 57.72000 40.58000 24.69000 28.26000 396.0 393.0
## 6 427.1236 423.4733 76.08643 66.38012 43.44486 25.77588 27.17419 375.0 372.5
## RFEB RFAB RFHD RFMLD RFAPD LTML LTPB LTMLD LTAPD RTML
## 1 85.0 75.85 48.36 28.94 33.01 353.0617 69.3506 21.46419 26.35661 365.5
## 2 73.0 62.00 44.09 24.11 24.41 283.0000 69.3506 19.12000 23.73000 281.5
## 3 77.5 67.35 44.17 27.18 28.39 330.0000 73.0000 21.72000 26.96000 331.5
## 4 77.5 66.95 43.77 25.18 31.50 312.0000 69.0000 21.76000 24.56000 320.5
## 5 70.0 57.88 41.13 24.79 26.88 306.0000 64.0000 19.44000 23.43000 305.0
## 6 68.0 58.90 40.80 23.76 24.40 298.0000 65.0000 17.75000 24.64000 293.5
## RTPB RTMLD RTAPD BIB LIBL RIBL LAcH RAcH
## 1 77.00000 22.41 26.75 274.0000 162.0000 161.00 52.56000 52.24000
## 2 69.34963 20.04 21.35 266.0000 151.0209 151.13 49.21000 50.18000
## 3 73.00000 21.07 25.36 262.3938 151.0209 151.13 48.88727 48.93736
## 4 71.50000 21.98 24.41 240.5000 151.0209 145.00 50.30000 49.80000
## 5 69.34963 20.30 21.45 285.5000 160.0000 155.00 45.60000 44.83000
## 6 62.50000 17.05 24.15 262.3938 151.0209 151.13 45.70000 47.04000sex_test## [1] 0 0 1 0 1 1 1 0 0 1
## Levels: 0 1# Feature Scaling: necesitamos escalar las variables en train y test para que esten entre -1 y 1
metricsgen_train = scale(metricsgen_train)
metricsgen_test = scale(metricsgen_test)
# Aplicamos el algoritmo SVM al dataset de train con la variable de prediccion sex_train
#install.packages('e1071')
library(e1071)
classifier = svm(formula = sex_train ~ .,
data = metricsgen_train,
type = 'C-classification',
kernel = 'linear')
# predecimos los resultados de test con el classifier creado
y_pred = predict(classifier, newdata = metricsgen_test)
# creamos la matriz de confusion para ver si ha funcionado
cm = table(sex_test, y_pred)
cm## y_pred
## sex_test 0 1
## 0 4 1
## 1 2 3Leyendo a la matriz de confusion, tenemos a 4 verdaderos positivos, 1 falso positivo, 2 falsos negativos y 3 verdaderos negativos. La accuracy es entonces de 7/10 o sea que en 7 casos de los 10, los resultados creados con SVM son los mismos que los que los autores han usado para predecir el genero de las observaciones donde el genero faltaba. Ahora el problema es que no sabiendo si las observaciones marcadas con ? en la columna Sex apartienen de verdad al genero asignado, es dificil saber si nuestro algoritmo a sido eficiente.
Apartado 4
El data set que tenemos contiene muchas variables que podemos cruzar en varios modos para hacernos unas preguntas.
- Unificacion de los formatos de la variable $Age:
Los catagorias de edad no tienen la misma composicion (20‐22; 22‐25; 25‐30; 30‐40; 40‐50; 50+) o sea que entre 20 y 30 anos tenemos a 3 categorias diferentes (20‐22; 22‐25; 25‐30) y todas las edades superiores a 50 estan en una categoria (50+). Mirando las observaciones por categoria, vemos que hay aun mas dispersion de edades:
library(ggplot2)
barplot(sort(table(data1$Age),decreasing = FALSE),horiz=TRUE, las=2, cex.names=.6)
Si si usa la variable Age, nos conviene unificar los valores de las observaciones, estableciendo unas categorias fijas, por ejempo <25 (dado que hay alcunas observaciones que estan en la categoria empezando con 18), 25-30, 30-40, 40-50, 50+
data2<-data1
data2$Age<-replace(data2$Age, data2$Age<25,"18-25")
data2$Age<-replace(data2$Age, data2$Age>=50,"50+")
data2$Age<-replace(data2$Age, data2$Age>=25 & data2$Age<30,"25-30")
data2$Age<-replace(data2$Age, data2$Age>=30 & data2$Age<40,"30-40")
data2$Age<-replace(data2$Age, data2$Age>=40 & data2$Age<50,"40-50")
#miramos la distribucion de las nuevas categorias de edades ahora
sort(table(data2$Age),decreasing = TRUE)##
## 30-40 40-50 25-30 50+ 18-25
## 576 443 250 157 112barplot(table(data2$Age), las=2, cex.names=.6)
Tenemos ahora a una distribucion similar a la normal con la mayoria de las observaciones aparteniendo a individuos entre 30 y 50 anos.
- Investigacion de la variable Location
Seria interesante también ver la distribucion de los orígenes de las observaciones (para un usuario no familiarizado con los sitios arqueologicos, en vez de usar la variable Nota que tiene el nombre del sitio arqueologico de origen, usaremos la variable Location)
sort(table(data2$Location),decreasing = TRUE)##
## Ohio, United States Germany
## 191 108
## Alaska, United States Japan
## 104 94
## Egypt Austria
## 81 78
## New Mexico, United States England, United Kingdom
## 73 72
## Kentucky, United States South Dakota, United States
## 61 61
## Aleutian Islands, Alaska, United States Sudan
## 56 53
## Italy Illinois, United States
## 51 42
## Belgium Utah, United States
## 41 35
## California, United States Philippine Islands
## 31 31
## Australia Peru
## 29 29
## Arizona, United States Argentina
## 23 21
## Chile New Jersey, United States
## 21 17
## Andaman Islands Madagascar
## 15 14
## Washington, United States Scotland, United Kingdom
## 14 13
## France Ecuador
## 12 11
## Canary Islands Colorado, United States
## 10 8
## Russia Democratic Republic of the Congo
## 8 6
## Arkansas, United States Solomon Islands
## 4 4
## South Africa Greenland
## 4 3
## China Hawaii, United States
## 2 2
## Indonesia Malaysia
## 1 1
## Papua New Guinea Tanzania
## 1 1
## Tasmania
## 1par(mar=c(4,10,4,4))
barplot(sort(table(data2$Location),decreasing = FALSE),horiz=TRUE, las=2, cex.names=.6)
length(unique(data2$Location))## [1] 45Vemos que la mayoria de las observaciones tienen origen en EEUU, seguidos por Egipto y Italia.
- Investigacion de la variable $Inst
Finalmente, nos gustaria ver un poco la distribucion de las observaciones por museo (variable Inst)
sort(table(data2$Inst),decreasing = TRUE)##
## NMNH AMNH CMNH MdH NM SfAP NHM MNdAE KU WOAC DC IRSN KSU
## 298 211 175 147 141 109 105 98 66 61 54 52 21par(mar=c(4,10,4,4))
barplot(sort(table(data2$Inst),decreasing = FALSE),horiz=TRUE, las=2, cex.names=.6)
Los museos donde se encuentran mas huesos observados son NMNH (Smithsonian National Museum of Natural History), AMNH (American Museum of Natural History), CMNH (Cleveland Museum of Natural History), MdH (Musee de l’Homme) y NM (Naturhistorisches Museum)
- Variable $Sex
sort(table(data2$Sex),decreasing = TRUE)##
## 0 1 0? 1?
## 985 543 5 5par(mar=c(4,10,4,4))
barplot(sort(table(data2$Sex),decreasing = FALSE),horiz=TRUE, las=2, cex.names=.6)
# porcentaje de observaciones donde no estan seguros del genero
10/1538## [1] 0.006501951Vemos que tenemos casi el doble de observaciones provenientes de individuos hombres que de las mujeres. El numero de observaciones donde los autores no estan seguros del genero es inferior a 1% (es de 0.65%). Hemos intentado en apartado 3 determinar el genero con un algoritmo de clasificacion SVM pero no llegamos a los mismos resultados que los autores (teniamos 30% de error)
- Frecuencias de las variables de tipo Elements (LHUM, RHUM, LRAD, RRAD, LFEM, RFEM, LTIB, RTIB, OSCX)
table(data2$LHUM)##
## 0 1
## 1383 155table(data2$RHUM)##
## 0 1
## 1413 125table(data2$LRAD)##
## 0 1
## 1322 216table(data2$RRAD)##
## 0 1
## 1343 195table(data2$LFEM)##
## 0 1
## 1430 108table(data2$RFEM)##
## 0 1
## 1439 99table(data2$LTIB)##
## 0 1
## 1405 133table(data2$RTIB)##
## 0 1
## 1405 133table(data2$OSCX)##
## 0 1
## 1507 31par(mfrow=c(4,2))
barplot(table(data2$LHUM),horiz=TRUE,xlab="LHUM")
barplot(table(data2$RHUM),horiz=TRUE,xlab="RHUM")
barplot(table(data2$LRAD),horiz=TRUE,xlab="LRAD")
barplot(table(data2$RRAD),horiz=TRUE,xlab="RRAD")
barplot(table(data2$LFEM),horiz=TRUE,xlab="LFEM")
barplot(table(data2$RFEM),horiz=TRUE,xlab="RFEM")
barplot(table(data2$LTIB),horiz=TRUE,xlab="LTIB")
barplot(table(data2$RTIB),horiz=TRUE,xlab="RTIB")
barplot(table(data2$OSCX),horiz=TRUE,xlab="OSCX")
Podemos ver que en todas las variables predomina la presencia del elemento y no su ausencia.
En cuanto a las variables del tipo Metrics, podemos ver su distribucion con unos Boxplot. Lo intuitivo puede ser mirar al lado las variables los huesos de izquierda y su equivalente de derecha
# ponemos las metricas en un dataset a parte
metrics <- data2[,c(18:60)]
# reordenamos la metrica para tener la izquierda al lado de la derecha para cada variable
metrics <- metrics[, c(1,6,2,7,3,8,4,9,5,10,11,14,12,15,13,16,17,24,18,25,19,26,20,27,21,28,22,29,23,30,31,35,32,36,33,37,34,38,40,41,42,43,39)]
boxplot(metrics,las=2, cex.names=.2)
Podemos ver ya que hay siempre una cierta disproporcion entre las variables de la izquierda y de la derecha, que sea de tamano promedio o de la varianza.
Apartado 5
Este data set permite múltiples cálculos probabilísticos, dada la gran cantidad de variables y distribuciones que contiene.
Si se selecciona un individuo al azar, ¿Cual es la probabilidad de que su achura bi-ilíaca (BIB) sea mayor que 270 mm?
#Primero se prepara un data set con los datos de la variable BIB (quitando los valores NA):
data_BIB_todaslasvariables <- data1[!is.na(data1$BIB),]
data_BIB <- data_BIB_todaslasvariables$BIB
#Se comprueba normalidad
ks.test(data_BIB, 'pnorm')## Warning in ks.test(data_BIB, "pnorm"): ties should not be present for the
## Kolmogorov-Smirnov test##
## One-sample Kolmogorov-Smirnov test
##
## data: data_BIB
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided#como hay normalidad, se calcula media y desviación típica de la distribución para calcular las probabilidades.
mean(data_BIB)## [1] 262.3958sd(data_BIB)## [1] 18.27456pnorm(q=270, mean = mean(data_BIB), sd = sd(data_BIB), lower.tail = FALSE)## [1] 0.3386662La probabilidad es del 33,87%.
Si se selecciona una mujer al azar, ¿Cual es la probabilidad de que su achura bi-ilíaca (BIB) sea mayor que 270 mm?
#Se prepara un nuevo data set, esta vez solo con mujeres, de la variable BIB.
BIB_mujeres <- data1[data1$Sex == '1', ]
data_BIB_variablesmujeres <- BIB_mujeres[!is.na(BIB_mujeres$BIB),]
data_BIB_mujeres <- data_BIB_variablesmujeres$BIB
#Se comprueba normalidad
ks.test(data_BIB_mujeres, 'pnorm')## Warning in ks.test(data_BIB_mujeres, "pnorm"): ties should not be present for
## the Kolmogorov-Smirnov test##
## One-sample Kolmogorov-Smirnov test
##
## data: data_BIB_mujeres
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided#como hay normalidad, se calcula media y desviación típica de la distribución para calcular las probabilidades.
mean(data_BIB_mujeres)## [1] 256.4375sd(data_BIB_mujeres)## [1] 17.81262pnorm(q=270, mean = mean(data_BIB_mujeres), sd = sd(data_BIB_mujeres), lower.tail = FALSE)## [1] 0.2232096La probabilidad es del 22,32%.
Si se selecciona un hombre al azar, ¿Cual es la probabilidad de que su achura bi-ilíaca (BIB) esté comprendida entre los 260 y los 280 mm?
#Se prepara otro data set, esta vez solo con hombres, de la variable BIB.
BIB_hombres <- data1[data1$Sex == '0', ]
data_BIB_variableshombres <- BIB_hombres[!is.na(BIB_hombres$BIB),]
data_BIB_hombres <- data_BIB_variableshombres$BIB
#Se comprueba normalidad
ks.test(data_BIB_hombres, 'pnorm')## Warning in ks.test(data_BIB_hombres, "pnorm"): ties should not be present for
## the Kolmogorov-Smirnov test##
## One-sample Kolmogorov-Smirnov test
##
## data: data_BIB_hombres
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided#como hay normalidad, se usa la media y desviación típica de la distribución para calcular la probabilidad de interés. Se emplea la fórmula siguiente:
1-pnorm(q=260, mean = mean(data_BIB_hombres), sd = sd(data_BIB_hombres), lower.tail = TRUE)-pnorm(q=280, mean = mean(data_BIB_hombres), sd = sd(data_BIB_hombres), lower.tail = FALSE)## [1] 0.4160365La probabilidad es del 41,60%.
Apartado 6
Antes de empezar, para poder realizar las regresiones lineales, es necesario modificar el data frame para excluir los valores NA. Para ello se itera en el data set, sustituyendo cada NA por la media de su columna:
#Se importa zoo para realizar la sustitución rápidamente
library(zoo)
metrics2 <- na.aggregate(metrics)
#Dado que el pool de variables es muy grande, el dataframe metrics2 creado anteriormente se divide en dos.
metrics2_1 <- metrics2[, c(1,6,2,7,3,8,4,9,5,10,11,14,12,15,13,16,17,24,18,25,19,26,20,27)]
metrics2_2 <- metrics2[, c(21,28,22,29,23,30,31,35,32,36,33,37,34,38,40,41,42,43,39)]
#Se observa si hay correlación entre variables.
library(corrplot)## corrplot 0.92 loadedcorrplot(cor(metrics2_1), method="number",order="hclust")
corrplot(cor(metrics2_2), method="number",order="hclust")
Ademés de la evidente correlación entre lados izquierdo y derecho para las mismas medidas, existe en general una alta correlación entre muchos pares. Como es dificil leer el resultado en una tabla con tantas variables, se hace de nuevo una separación en tablas mas pequeñas.
metrics2_1_1 <- metrics2_1[, c(1,2,3,4,5,6,7,8)]
metrics2_1_2 <- metrics2_1[, c(9,10,11,12,13,14,15,16)]
metrics2_1_3 <- metrics2_1[, c(17,18,19,20,21,22,23,24)]
metrics2_2 <- metrics2[, c(21,28,22,29,23,30,31,35,32,36,33,37,34,38,40,41,42,43,39)]
metrics2_2_1 <- metrics2_2[, c(1,2,3,4,5,6,7,8)]
metrics2_2_2 <- metrics2_2[, c(9,10,11,12,13,14,15,16,17,18,19)]
corrplot(cor(metrics2_1_1), method="number",order="hclust")
corrplot(cor(metrics2_1_2), method="number",order="hclust")
corrplot(cor(metrics2_1_3), method="number",order="hclust")
corrplot(cor(metrics2_2_1), method="number",order="hclust")
corrplot(cor(metrics2_2_2), method="number",order="hclust")
La correlacion mas alta (entre variables que no sean gemelas), se da entre los pares RHEB-RHHD, RFML-RFBL, LFML-LFBL y LFAB-RFEB.
Ahora se realizan las pruebas de regresión lineal entre estos pares:
Regresion1<-lm(RHEB~RHHD, data=metrics2)
summary(Regresion1)##
## Call:
## lm(formula = RHEB ~ RHHD, data = metrics2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.260 -1.766 0.000 1.694 10.953
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.8584 0.7948 14.92 <2e-16 ***
## RHHD 1.0775 0.0184 58.56 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.926 on 1536 degrees of freedom
## Multiple R-squared: 0.6906, Adjusted R-squared: 0.6904
## F-statistic: 3429 on 1 and 1536 DF, p-value: < 2.2e-16plot(RHEB~RHHD, data=metrics2)
abline(11.86, 1.08)
abline(Regresion1)
Regresion2<-lm(RFML~RFBL, data=metrics2)
summary(Regresion2)##
## Call:
## lm(formula = RFML ~ RFBL, data = metrics2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -94.007 -1.610 -0.142 1.146 97.314
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.388530 1.747054 5.946 3.39e-09 ***
## RFBL 0.984936 0.004133 238.323 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.963 on 1536 degrees of freedom
## Multiple R-squared: 0.9737, Adjusted R-squared: 0.9737
## F-statistic: 5.68e+04 on 1 and 1536 DF, p-value: < 2.2e-16plot(RFML~RFBL, data=metrics2)
abline(10.39, 0.94)
abline(Regresion2)
Regresion3<-lm(LFML~LFBL, data=metrics2)
summary(Regresion3)##
## Call:
## lm(formula = LFML ~ LFBL, data = metrics2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -101.925 -1.327 0.000 1.193 79.624
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.203046 1.506410 5.445 6.01e-08 ***
## LFBL 0.989251 0.003548 278.792 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.218 on 1536 degrees of freedom
## Multiple R-squared: 0.9806, Adjusted R-squared: 0.9806
## F-statistic: 7.773e+04 on 1 and 1536 DF, p-value: < 2.2e-16plot(LFML~LFBL, data=metrics2)
abline(8.20, 0.99)
abline(Regresion3)
Regresion4<-lm(LFAB~RFEB, data=metrics2)
summary(Regresion4)##
## Call:
## lm(formula = LFAB ~ RFEB, data = metrics2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.026 -1.499 0.000 1.526 16.209
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.16829 0.94165 5.489 4.73e-08 ***
## RFEB 0.80243 0.01231 65.180 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.858 on 1536 degrees of freedom
## Multiple R-squared: 0.7345, Adjusted R-squared: 0.7343
## F-statistic: 4248 on 1 and 1536 DF, p-value: < 2.2e-16plot(LFAB~RFEB, data=metrics2)
abline(5.17, 0.80)
abline(Regresion4)
Para la regresión múltiple, hemos decidido tratar de predecir la anchura epicondilar del húmero derecho (RHEB) a partir de las otras 7 variables que también han sido usadas en la regresión simple. Se observa el modelo y también se usa el método stepwise para comprobar cuales de las variables son las mejores predictoras.
mrm_metrics2 <- lm(metrics2$RHEB~metrics2$RHHD+metrics2$RFML+metrics2$RFBL+metrics2$LFML+metrics2$LFBL+metrics2$LFAB+metrics2$RFEB)
summary(mrm_metrics2)##
## Call:
## lm(formula = metrics2$RHEB ~ metrics2$RHHD + metrics2$RFML +
## metrics2$RFBL + metrics2$LFML + metrics2$LFBL + metrics2$LFAB +
## metrics2$RFEB)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.7763 -1.7429 -0.0294 1.7024 9.6017
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.024e+00 1.055e+00 2.868 0.00419 **
## metrics2$RHHD 7.098e-01 3.199e-02 22.186 < 2e-16 ***
## metrics2$RFML 1.490e-03 1.538e-02 0.097 0.92281
## metrics2$RFBL -1.637e-06 1.421e-02 0.000 0.99991
## metrics2$LFML -2.025e-02 1.762e-02 -1.150 0.25051
## metrics2$LFBL 3.018e-02 1.689e-02 1.787 0.07407 .
## metrics2$LFAB 2.028e-02 2.808e-02 0.722 0.47029
## metrics2$RFEB 2.431e-01 2.801e-02 8.678 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.754 on 1530 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7258
## F-statistic: 582.1 on 7 and 1530 DF, p-value: < 2.2e-16step(object=mrm_metrics2,direction ="both", trace=1)## Start: AIC=3124.18
## metrics2$RHEB ~ metrics2$RHHD + metrics2$RFML + metrics2$RFBL +
## metrics2$LFML + metrics2$LFBL + metrics2$LFAB + metrics2$RFEB
##
## Df Sum of Sq RSS AIC
## - metrics2$RFBL 1 0.0 11605 3122.2
## - metrics2$RFML 1 0.1 11605 3122.2
## - metrics2$LFAB 1 4.0 11609 3122.7
## - metrics2$LFML 1 10.0 11615 3123.5
## <none> 11605 3124.2
## - metrics2$LFBL 1 24.2 11629 3125.4
## - metrics2$RFEB 1 571.2 12176 3196.1
## - metrics2$RHHD 1 3733.5 15338 3551.2
##
## Step: AIC=3122.18
## metrics2$RHEB ~ metrics2$RHHD + metrics2$RFML + metrics2$LFML +
## metrics2$LFBL + metrics2$LFAB + metrics2$RFEB
##
## Df Sum of Sq RSS AIC
## - metrics2$RFML 1 0.4 11605 3120.2
## - metrics2$LFAB 1 4.0 11609 3120.7
## - metrics2$LFML 1 10.1 11615 3121.5
## <none> 11605 3122.2
## - metrics2$LFBL 1 24.4 11629 3123.4
## + metrics2$RFBL 1 0.0 11605 3124.2
## - metrics2$RFEB 1 571.6 12176 3194.1
## - metrics2$RHHD 1 3733.5 15338 3549.2
##
## Step: AIC=3120.23
## metrics2$RHEB ~ metrics2$RHHD + metrics2$LFML + metrics2$LFBL +
## metrics2$LFAB + metrics2$RFEB
##
## Df Sum of Sq RSS AIC
## - metrics2$LFAB 1 3.6 11609 3118.7
## - metrics2$LFML 1 9.9 11615 3119.5
## <none> 11605 3120.2
## - metrics2$LFBL 1 24.6 11630 3121.5
## + metrics2$RFML 1 0.4 11605 3122.2
## + metrics2$RFBL 1 0.3 11605 3122.2
## - metrics2$RFEB 1 686.5 12292 3206.6
## - metrics2$RHHD 1 3766.0 15371 3550.5
##
## Step: AIC=3118.71
## metrics2$RHEB ~ metrics2$RHHD + metrics2$LFML + metrics2$LFBL +
## metrics2$RFEB
##
## Df Sum of Sq RSS AIC
## - metrics2$LFML 1 10.6 11619 3118.1
## <none> 11609 3118.7
## + metrics2$LFAB 1 3.6 11605 3120.2
## - metrics2$LFBL 1 26.7 11635 3120.2
## + metrics2$RFBL 1 0.0 11609 3120.7
## + metrics2$RFML 1 0.0 11609 3120.7
## - metrics2$RFEB 1 1074.8 12683 3252.9
## - metrics2$RHHD 1 4083.6 15692 3580.3
##
## Step: AIC=3118.11
## metrics2$RHEB ~ metrics2$RHHD + metrics2$LFBL + metrics2$RFEB
##
## Df Sum of Sq RSS AIC
## <none> 11619 3118.1
## + metrics2$LFML 1 10.6 11609 3118.7
## + metrics2$LFAB 1 4.3 11615 3119.5
## + metrics2$RFML 1 1.2 11618 3120.0
## + metrics2$RFBL 1 0.8 11618 3120.0
## - metrics2$LFBL 1 102.7 11722 3129.6
## - metrics2$RFEB 1 1072.5 12692 3251.9
## - metrics2$RHHD 1 4078.2 15697 3578.8##
## Call:
## lm(formula = metrics2$RHEB ~ metrics2$RHHD + metrics2$LFBL +
## metrics2$RFEB)
##
## Coefficients:
## (Intercept) metrics2$RHHD metrics2$LFBL metrics2$RFEB
## 2.9041 0.7153 0.0120 0.2550Según el método stepwise, las mejores variables de este conjunto para predecir RHEB son RHHD, LFBL y RFEB.
Apartado 7
Está documentado en Antropología Física desde hace mucho tiempo el hecho de que existen diferencias significativas entre sexos para una gran cantidad de variables óseas lineales. Generalmente, el tamaño de los huesos masculinos es mayor que el de los femeninos. En el data set que hemos empleado hay multitud de ejemplos de este fenómeno, pues suele darse abundantemente en variables de huesos largos.
Por esta razón hemos decidido usar una variable de la cadera (formada por huesos planos en vez de largos), y asi dilucidar si el fenómeno de dimorfismo sexual ocurre también en esta region y en esta población. La variable elegida es la Anchura Bi-ilíaca (BIB).
El primer paso ha sido quitar los casos del data set que tienen los valores “0?” y “1?” para la variable Sex (es decir, los casos cuyo sexo no es seguro).
data_ANOVA<-data1[!(data1$Sex=="0?" | data1$Sex=="1?"),]Después se ha cambiado la variable Sex de tipo “caracter” a tipo “factor” (necesario para la ANOVA):
data_ANOVA$Sex <- as.factor(data_ANOVA$Sex)Ahora se puede comenzar con la prueba ANOVA. Hay que recordar que para que la anova pueda implementarse, primero hay que demostrar normalidad entre las dos poblaciones que se comparen (en este caso, hombres y mujeres), asi como igualdad de varianzas en sus distribuciones.
ANOVA
Se crea un data frame nuevo que contenga solo las variables Sex y BIB, a partir del data frame preparado en los pasos anteriores, y que además no contenga casos con NA en dichas variables:
ANOVA_1 <- data.frame(data_ANOVA$Sex, data_ANOVA$BIB)
ANOVA_1 <- na.omit(ANOVA_1)Se prueba normalidad para ambos sexos en la variable BIB:
library(nortest)
by(data=ANOVA_1,INDICES = ANOVA_1$data_ANOVA.Sex,FUN=function(x){lillie.test(x$data_ANOVA.BIB)})## ANOVA_1$data_ANOVA.Sex: 0
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$data_ANOVA.BIB
## D = 0.038449, p-value = 0.002299
##
## ------------------------------------------------------------
## ANOVA_1$data_ANOVA.Sex: 1
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$data_ANOVA.BIB
## D = 0.046259, p-value = 0.00989La prueba indica que hay normalidad en ambos grupos. Ahora se hace la prueba de igualdad de varianzas:
fligner.test(ANOVA_1$data_ANOVA.BIB~ANOVA_1$data_ANOVA.Sex,ANOVA_1)##
## Fligner-Killeen test of homogeneity of variances
##
## data: ANOVA_1$data_ANOVA.BIB by ANOVA_1$data_ANOVA.Sex
## Fligner-Killeen:med chi-squared = 0.043118, df = 1, p-value = 0.8355La prueba indica que hay igualdad de varianzas, por ello, se continua con la ANOVA
anova<-aov(ANOVA_1$data_ANOVA.BIB~ANOVA_1$data_ANOVA.Sex,data=ANOVA_1)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## ANOVA_1$data_ANOVA.Sex 1 28664 28664 90.87 <2e-16 ***
## Residuals 1457 459601 315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El resultado de la ANOVA parece confirmar que existen diferencias significativas en la anchura bi-ilíaca entre hombres y mujeres de esta colección.
Apartado 8
Estos datos han permitido obtener una gran cantidad de resultados y llevar a cabo muchas pruebas de interés gracias a, principalmente, tres factores del data set. El primero es la abundancia en si del registro: la gran cantidad de casos distintos que recoge y la gran variedad de variables tomadas. Esto ha permitido seleccionar las variables mas indicadas para cada prueba que se ha realizado, teniendo siempre disponibles los suficientes individuos para llevarla a cabo. En segundo lugar, las variables que recoge son de distinto tipo, no teniendo únicamente variables continuas como cabe esperar de un data set osteológico, sino tambien de tipo factor (sexo, yacimiento) o numéricas discretas (edad). Esto permite realizar muchas pruebas de alto valor estadístico, como la ANOVA. Por último, dada la naturaleza de los datos que recoge, este data set tiene muchas variables “por parejas”, es decir, que de cada medida ha recogido el valor tanto del lado derecho como del izquierdo. Esto permite realizar estudios de lateralidad, sustituir una variable por su homóloga hermana si no hubiese suficientes casos o se hubiese perdido información, o , si se supiese su lado preferencial (diestro o zurdo) hacer análisis de la influecia de esta variable sobre el desarrollo diferencial de cada lado del esqueleto.
Hemos podido obtener algunas conclusiones intermediarias con el estudio del dataset que hemos hecho. No obstante ello, la cantidad de variables y observaciones presenta muchas mas oportunidades par un estudio mas produndizado que dependeria de la direccion que se quiere tomar (unos huesos especificos? de una origen, genero, tamano especificos? haciendo unos algoritmos predictivos (como en el caso del SVM), usar otro algoritmo (por ejemplo regresion logistica o KNN)? en vez de tener un dataset de train y uno de test, usar también un subset del train para validacion?)